2020-07-04 - Environment Change Rate Sweeps - Experiment Analysis
Analysis Dependencies
library(ggplot2) # (Wickham, 2016)
library(tidyr) # (Wickham and Henry, 2020)
library(dplyr) # (Wickham et al., 2020)
library(reshape2) # (Wickham, 2007)
library(cowplot) # (Wilke, 2019)
library(hexbin) # (Carr, Lewin-Kog, Maechler, and Sarkar, 2020)We conducted these analyses using the following computing environment:
print(version)## _
## platform x86_64-apple-darwin15.6.0
## arch x86_64
## os darwin15.6.0
## system x86_64, darwin15.6.0
## status
## major 3
## minor 6.2
## year 2019
## month 12
## day 12
## svn rev 77560
## language R
## version.string R version 3.6.2 (2019-12-12)
## nickname Dark and Stormy NightSetup
theme_set(theme_cowplot())
alpha <- 0.05
data_path <- "./data/agg_data.csv"
agg_data <- read.csv(data_path, na.strings="NONE")
agg_data$BIT_FLIP_PROB <- as.factor(agg_data$BIT_FLIP_PROB)
agg_data$DRIFT <- agg_data$TOURNAMENT_SIZE==1
# Label
chg_rate_label <- function(mag, interval, drift) {
if (drift) { return("drift") }
else if (interval == 0) { return("0") }
else { return(paste(mag, interval, sep="/")) }
}
agg_data$chg_rate_label <- factor(mapply(chg_rate_label,
agg_data$CHANGE_MAGNITUDE,
agg_data$CHANGE_FREQUENCY,
agg_data$DRIFT),
levels=c("drift", "0", "1/256", "1/128",
"1/64", "1/32", "1/16", "1/8",
"1/4", "1/2", "1/1", "2/1",
"4/1", "8/1", "16/1", "32/1",
"64/1", "128/1", "256/1",
"512/1", "1024/1", "2048/1",
"4096/1"))
# Divide up the data into convenient partitions
data_gradient_phase0 <-
filter(agg_data,
GRADIENT_MODEL==1 & evo_phase == 0 & update == 50000)
data_gradient_phase1 <-
filter(agg_data,
GRADIENT_MODEL==1 & evo_phase == 1 & update == 60000)
data_nk_phase0 <-
filter(agg_data,
GRADIENT_MODEL==0 & evo_phase == 0 & update == 50000)
data_nk_phase1 <-
filter(agg_data,
GRADIENT_MODEL==0 & evo_phase == 1 & update == 60000)Evolution Experiment
We evolved populations for 50,000 generations under a range of environmental change rates.
Fitnesses
Note issue with fitnesses for changing environments.
Gradient Fitness Model
ggplot(data_gradient_phase0,
aes(x=chg_rate_label, y=fitness, color=chg_rate_label)) +
geom_boxplot() +
xlab("Environmental Change") +
scale_y_continuous(name="Fitness (of best organism)") +
ggtitle("Gradient Fitness Model") +
theme(legend.position="none")
NK Fitness Model
ggplot(data_nk_phase0,
aes(x=chg_rate_label, y=fitness, color=chg_rate_label)) +
geom_boxplot() +
xlab("Env. bits flipped per generation") +
scale_y_continuous(name="Fitness (of best organism)") +
ggtitle("NK Fitness Model") +
theme(legend.position="none")
Coding Sites
Gradient Fitness Model
ggplot(data_gradient_phase0,
aes(x=chg_rate_label, y=coding_sites, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(aes(color=chg_rate_label), alpha=0.1) +
xlab("Environmental Change") +
scale_y_continuous(name="Coding Sites (in best genotype)",
limits=c(0, 130),
breaks=seq(0, 130, 20)) +
ggtitle("Gradient Fitness Model") +
theme(legend.position="none")
kt <- kruskal.test(formula = coding_sites ~ chg_rate_label, data=data_gradient_phase0)
kt##
## Kruskal-Wallis rank sum test
##
## data: coding_sites by chg_rate_label
## Kruskal-Wallis chi-squared = 811.32, df = 10, p-value < 2.2e-16if (kt$p.value <= alpha) {
wt <- pairwise.wilcox.test(x=data_gradient_phase0$coding_sites,
g=data_gradient_phase0$chg_rate_label,
exact=FALSE,
p.adjust.method = "bonferroni")
wt
}##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: data_gradient_phase0$coding_sites and data_gradient_phase0$chg_rate_label
##
## drift 0 1/256 1/128 1/64 1/32 1/16 1/8 1/4
## 0 1.0000 - - - - - - - -
## 1/256 < 2e-16 < 2e-16 - - - - - - -
## 1/128 < 2e-16 < 2e-16 3.2e-07 - - - - - -
## 1/64 < 2e-16 < 2e-16 < 2e-16 4.6e-13 - - - - -
## 1/32 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.0194 - - - -
## 1/16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 6.4e-10 0.0178 - - -
## 1/8 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.5e-14 2.5e-06 1.0000 - -
## 1/4 < 2e-16 < 2e-16 < 2e-16 < 2e-16 2.4e-10 0.0020 1.0000 1.0000 -
## 1/2 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.0929 1.0000 0.0210 1.2e-05 0.0028
## 1/1 < 2e-16 < 2e-16 < 2e-16 5.1e-06 0.3774 6.2e-07 1.6e-14 < 2e-16 1.8e-14
## 1/2
## 0 -
## 1/256 -
## 1/128 -
## 1/64 -
## 1/32 -
## 1/16 -
## 1/8 -
## 1/4 -
## 1/2 -
## 1/1 4.4e-06
##
## P value adjustment method: bonferroniNK Fitness Model
ggplot(data_nk_phase0,
aes(x=chg_rate_label, y=coding_sites, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
scale_y_continuous(name="# coding sites in best organism",
limits=c(0, 130),
breaks=seq(0, 130, 20)) +
ggtitle("NK Fitness Model") +
theme(legend.position="none")
kt <- kruskal.test(formula = coding_sites ~ chg_rate_label, data=data_nk_phase0)
kt##
## Kruskal-Wallis rank sum test
##
## data: coding_sites by chg_rate_label
## Kruskal-Wallis chi-squared = 536.31, df = 10, p-value < 2.2e-16if (kt$p.value <= alpha) {
wt <- pairwise.wilcox.test(x=data_nk_phase0$coding_sites,
g=data_nk_phase0$chg_rate_label,
exact=FALSE,
p.adjust.method = "bonferroni")
wt
}##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: data_nk_phase0$coding_sites and data_nk_phase0$chg_rate_label
##
## drift 0 1/1 2/1 4/1 8/1 16/1 32/1 64/1
## 0 1.00000 - - - - - - - -
## 1/1 < 2e-16 < 2e-16 - - - - - - -
## 2/1 < 2e-16 < 2e-16 1.00000 - - - - - -
## 4/1 < 2e-16 < 2e-16 0.35262 1.00000 - - - - -
## 8/1 < 2e-16 < 2e-16 0.00048 0.00592 1.00000 - - - -
## 16/1 < 2e-16 < 2e-16 6.4e-07 1.5e-05 0.23057 1.00000 - - -
## 32/1 < 2e-16 < 2e-16 5.4e-11 1.7e-09 6.0e-05 0.00551 1.00000 - -
## 64/1 < 2e-16 < 2e-16 1.3e-13 5.8e-12 1.0e-06 0.00010 0.11885 1.00000 -
## 128/1 < 2e-16 < 2e-16 9.8e-14 7.3e-12 9.5e-07 7.3e-05 0.10450 1.00000 1.00000
## 256/1 < 2e-16 < 2e-16 0.00211 0.03518 1.00000 1.00000 1.00000 0.01884 0.00082
## 128/1
## 0 -
## 1/1 -
## 2/1 -
## 4/1 -
## 8/1 -
## 16/1 -
## 32/1 -
## 64/1 -
## 128/1 -
## 256/1 0.00094
##
## P value adjustment method: bonferroniGenome Length
Gradient Fitness Model
ggplot(data_gradient_phase0,
aes(x=chg_rate_label, y=genome_length, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
ylab("Genome Length") +
ylim(0, 1025) +
ggtitle("Gradient Fitness Model") +
theme(legend.position="none")
NK Fitness Model
ggplot(data_nk_phase0,
aes(x=chg_rate_label, y=genome_length, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
ylab("Genome Length") +
ylim(0, 1024) +
ggtitle("NK Fitness Model") +
theme(legend.position="none")
Neutral Sites
Gradient Fitness Model
ggplot(data_gradient_phase0,
aes(x=chg_rate_label, y=neutral_sites, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
ylab("Neutral Sites") +
ylim(-1, 1025) +
ggtitle("Gradient Fitness Model") +
theme(legend.position="none")
NK Fitness Model
ggplot(data_nk_phase0,
aes(x=chg_rate_label, y=neutral_sites, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
ylab("Neutral Sites") +
ylim(-1, 1025) + # jitter can jitter below 0
ggtitle("NK Fitness Model") +
theme(legend.position="none")
Transfer experiment
We lock-down genetic architectures and evolve for 10,000 generations in a random static environment.
Fitnesses
Gradient Fitness Model
ggplot(data_gradient_phase1,
aes(x=chg_rate_label, y=fitness, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
scale_y_continuous(name="Fitness") +
ggtitle("Gradient Fitness Model (transfer)") +
theme(legend.position="none")
kt <- kruskal.test(formula = fitness ~ chg_rate_label, data=data_gradient_phase1)
kt##
## Kruskal-Wallis rank sum test
##
## data: fitness by chg_rate_label
## Kruskal-Wallis chi-squared = 742.4, df = 10, p-value < 2.2e-16if (kt$p.value <= alpha) {
wt <- pairwise.wilcox.test(x=data_gradient_phase1$fitness,
g=data_gradient_phase1$chg_rate_label,
exact=FALSE,
p.adjust.method = "bonferroni")
wt
}##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: data_gradient_phase1$fitness and data_gradient_phase1$chg_rate_label
##
## drift 0 1/256 1/128 1/64 1/32 1/16 1/8 1/4
## 0 1.00000 - - - - - - - -
## 1/256 < 2e-16 < 2e-16 - - - - - - -
## 1/128 < 2e-16 < 2e-16 2.0e-05 - - - - - -
## 1/64 < 2e-16 < 2e-16 < 2e-16 1.1e-09 - - - - -
## 1/32 < 2e-16 < 2e-16 < 2e-16 4.8e-15 1.00000 - - - -
## 1/16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.1e-06 0.00526 - - -
## 1/8 < 2e-16 < 2e-16 < 2e-16 < 2e-16 2.0e-09 2.4e-05 1.00000 - -
## 1/4 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.00298 1.00000 1.00000 1.00000 -
## 1/2 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.17302 1.00000 0.12693 0.00076 1.00000
## 1/1 < 2e-16 < 2e-16 < 2e-16 7.6e-05 1.00000 0.00153 7.7e-11 1.4e-13 8.6e-07
## 1/2
## 0 -
## 1/256 -
## 1/128 -
## 1/64 -
## 1/32 -
## 1/16 -
## 1/8 -
## 1/4 -
## 1/2 -
## 1/1 3.0e-05
##
## P value adjustment method: bonferroniNK Fitness Model
ggplot(data_nk_phase1,
aes(x=chg_rate_label, y=fitness, color=chg_rate_label)) +
geom_boxplot() +
geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
scale_y_continuous(name="Fitness") +
ggtitle("NK Fitness Model (transfer)") +
theme(legend.position="none")
kt <- kruskal.test(formula = fitness ~ chg_rate_label, data=data_nk_phase1)
kt##
## Kruskal-Wallis rank sum test
##
## data: fitness by chg_rate_label
## Kruskal-Wallis chi-squared = 449.57, df = 10, p-value < 2.2e-16if (kt$p.value <= alpha) {
wt <- pairwise.wilcox.test(x=data_nk_phase1$fitness,
g=data_nk_phase1$chg_rate_label,
exact=FALSE,
p.adjust.method = "bonferroni")
wt
}##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: data_nk_phase1$fitness and data_nk_phase1$chg_rate_label
##
## drift 0 1/1 2/1 4/1 8/1 16/1 32/1 64/1 128/1
## 0 1.0000 - - - - - - - - -
## 1/1 2.8e-16 < 2e-16 - - - - - - - -
## 2/1 5.1e-14 < 2e-16 1.0000 - - - - - - -
## 4/1 < 2e-16 < 2e-16 1.0000 0.0701 - - - - - -
## 8/1 < 2e-16 < 2e-16 1.0000 0.0040 1.0000 - - - - -
## 16/1 < 2e-16 < 2e-16 1.0000 0.0055 1.0000 1.0000 - - - -
## 32/1 < 2e-16 < 2e-16 0.0073 3.4e-06 0.5241 1.0000 1.0000 - - -
## 64/1 < 2e-16 < 2e-16 3.5e-05 1.6e-10 0.0075 0.0210 0.4969 1.0000 - -
## 128/1 < 2e-16 < 2e-16 8.1e-06 4.2e-10 0.0019 0.0056 0.1172 1.0000 1.0000 -
## 256/1 < 2e-16 < 2e-16 1.0000 0.0141 1.0000 1.0000 1.0000 1.0000 0.0884 0.0132
##
## P value adjustment method: bonferroniFitness vs Number of Coding Sites
Gradient Fitness Model
ggplot(data_gradient_phase1,
aes(x=coding_sites, y=fitness)) +
geom_point(aes(color=chg_rate_label)) +
xlab("Coding Sites") +
ylab("Fitness") +
theme(legend.position="top") +
ggtitle("Gradient Fitness Model")
cor.test(x=data_gradient_phase1$coding_sites,
y=data_gradient_phase1$fitness,
method="spearman",
exact=FALSE)##
## Spearman's rank correlation rho
##
## data: data_gradient_phase1$coding_sites and data_gradient_phase1$fitness
## S = 8352633, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.9623472NK Fitness Model
ggplot(data_nk_phase1,
aes(x=coding_sites, y=fitness)) +
geom_point(aes(color=chg_rate_label)) +
xlab("Coding Sites") +
ylab("Fitness") +
theme(legend.position="top") +
ggtitle("NK Fitness Model")
cor.test(x=data_nk_phase1$coding_sites,
y=data_nk_phase1$fitness,
method="spearman",
exact=FALSE)##
## Spearman's rank correlation rho
##
## data: data_nk_phase1$coding_sites and data_nk_phase1$fitness
## S = 38428070, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8267704Publication Figures
Changing Environments Promote Modularity
ggplot(data_gradient_phase0,
aes(x=chg_rate_label, y=coding_sites)) +
geom_boxplot() +
# geom_jitter(alpha=0.1) +
xlab("Environmental Change") +
scale_y_continuous(name="Coding Sites\n(best genotype)",
limits=c(0, 130),
breaks=seq(0, 130, 20)) +
theme(legend.position="none") +
ggsave("./imgs/coding-sites-v-chg-rate-gradient-phase0.pdf", width=10, height=7)
Modularity Promotes Evolvability
ggplot(data_gradient_phase1,
aes(x=chg_rate_label, y=fitness)) +
geom_boxplot() +
xlab("Environmental Change") +
scale_y_continuous(name="Fitness") +
# ggtitle("Gradient Fitness Model (phase 2)") +
theme(legend.position="none") +
ggsave("./imgs/fitness-v-chg-rate-gradient-phase1.pdf", width=10, height=7)
ggplot(data_gradient_phase1,
aes(x=coding_sites, y=fitness)) +
geom_hex(alpha=0.95) +
scale_fill_continuous(type = "viridis",
trans="log",
name="Count",
breaks=c(1, 10, 100)) +
scale_x_continuous(name="Coding Sites",
limits=c(0, 130),
breaks=seq(0, 130, 20)) +
ylab("Fitness") +
ggsave("./imgs/coding-sites-v-fitness-gradient-phase1-hex.pdf", width=10, height=7)
ggplot(data_gradient_phase1,
aes(x=coding_sites, y=fitness)) +
geom_density_2d() +
stat_density_2d(aes(fill = ..level..),
geom = "polygon",
colour="white") +
geom_jitter(alpha=0.05) +
scale_fill_continuous(type = "viridis",
trans="log",
name="Count") +
ylab("Fitness") +
theme(legend.position = "none")
ggplot(data_gradient_phase1,
aes(x=coding_sites, y=fitness)) +
geom_point(alpha=0.95) +
scale_x_continuous(name="Coding Sites",
limits=c(0, 130),
breaks=seq(0, 130, 20)) +
ylab("Fitness")